function S = Smatdatenums(gt,h,wt)
% Difference with Smatdates: h already contains datenums for start and
% ending dates

t = size(gt,1);
nq = size(gt,2);

st = h(:,2); 
en = h(:,3); 
S0 =zeros(nq,nq); 
v = zeros(nq,1); 

strtdatenum = st;
enddatenum = en;
% if st(1)>19800101
%     Tyr  = floor(st/10000);
%     Tmth = floor(st/100) - Tyr*100;
%     Tday = mod(st,100);
%     strtdatenum = datenum(Tyr, Tmth, Tday);
%     Tyr  = floor(en/10000);
%     Tmth = floor(en/100) - Tyr*100;
%     Tday = mod(en,100);
%     enddatenum = datenum(Tyr, Tmth, Tday);
% else
%     Tyr  = floor(st/100);
%     Tmth = st - Tyr*100;
%     strtdatenum = datenum(Tyr, Tmth, 1);
%     Tyr  = floor(en/100);
%     Tmth = en - Tyr*100;
%     enddatenum = datenum(Tyr, Tmth, 1);
% end    

%Example: for four obs of which two perfectly overlap and are perfectly
%correlated, the calculation below produces 
%S = 2Var(r_i), i.e. precisely what's needed to take into account that only
%two effective observations: s.e. of mean is (1/2)Var(r_i) rather than
%(1/4)Var(r_i) 


for i = 1:t 
    
    %also tried a truncated window estimator with weights of one if there is
    %any overlap, but this one can produce non PSD S matrix
    %ovlp = find(en >= st(i) & st <= en(i));  %find overlapping return periods
    %S = S + gt(ovlp,:)'*repmat(gt(i,:),size(gt(ovlp,:),1),1);
    %%Test: with this it reproduces Cov(gt): S = S + gt(i,:)'*repmat(gt(i,:),size(gt(i,:),1),1);
    
    
%     %this corresponds to Bartlett estimator in Conley's Spatial GMM
%     mins = min(st(i),st); 
%     maxs = max(st(i),st); 
%     mine = min(en(i),en); 
%     maxe = max(en(i),en); 

    %this corresponds to Bartlett estimator in Conley's Spatial GMM
    mins = min(strtdatenum(i),strtdatenum); 
    maxs = max(strtdatenum(i),strtdatenum); 
    mine = min(enddatenum(i),enddatenum); 
    maxe = max(enddatenum(i),enddatenum); 
   
    
    d = 1-(mine-maxs)./(maxe-mins);  %equals zero for complete overlap, one for none, greater than one for interval inbetween
    
    w = 1-min(d,wt)/(wt); 
    
    %S = S + (repmat(w,1,nq).*gt)'*repmat(gt(i,:),t,1);
    
    %this specification constructs sum of weighted cross products of
    %gt(i,k) and gt(j,k) for each asset class k, but the correlations
    %between asset classes are taken only from sum of gt(i,:)'*gt(i,:)
    %i.e., no serial (including cross-serial correlation of the underlying
    %returns process, so that 
    % Cov(x_t + x_{t+1}, y_t + y_{t+1})/(2 \sigma_x \sigma_y) 
    % = Cov(x_t , y_t )/(\sigma_x \sigma_y) 
   
    v = v + max(0,diag( (repmat(w,1,nq).*gt)'*repmat(gt(i,:),t,1) ));   % the max(0,.) is because for large t you sometimes numerically get a tiny negative number
    
    S0 = S0 + gt(i,:)'*gt(i,:);
    
end    

v0 = diag(S0); 
sc = (v./v0).^(1/2); 
S = (1/t)*diag(sc)*S0*diag(sc); 

%S = S/t; 


